#### REPLICATION MATERIALS FOR ''BENCHMARKING PANDEMIC RESPONSE: HOW THE UK'S COVID-19 VACCINE ROLL-OUT AFFECTED POPULAR SUPPORT FOR THE EU''
###### Authors: Irene Rodríguez, Toni Rodon, Asli Unan, Lisa Herbig, Heike Klüver & Theresa Kuhn
###### Journal: British Journal of Political Science

###### FIGURES AND TABLES IN BODY OF MANUSCRIPT

# Loading libraries
if (!require("pacman")) install.packages("pacman")
pacman::p_load(kableExtra, marginaleffects, modelsummary, knitr, estimatr, tidyverse, cobalt, 
               knitr, RColorBrewer, janitor, lubridate, readstata13, countrycode, rdrobust, 
               magick, jtools, ggiraphExtra, sjPlot, gridExtra, data.table, pastecs, MatchIt, readxl)

# Set working directory
setwd("~/Replication materials")

load("Data/data.RData")

#### Figure 1: Weekly proportion of newspaper articles' titles about the UK and the EU that mention vaccination ####

euvsuk_vax <- read_csv("Data/ukvseu_vax.csv") # importing python scraping files

euvsuk_vax <- euvsuk_vax |> # manipulating table to fit visualisation needs
  pivot_wider(names_from = "topic",
              values_from = "attention")

# assigning colors 
colors <- c("European Union" = "blue", "United Kingdom" = "red")

# building plot for Figure 1
p_vax_attention <- ggplot(data = euvsuk_vax,
                          mapping = aes(x = date)) +
  geom_point(aes(y = `United Kingdom`, color = 'United Kingdom'), size = 1) +
  geom_smooth(aes(y = `United Kingdom`), color = 'red', method = "loess") +
  geom_point(aes(y = `European Union`, color = 'European Union'), size = 1) +
  geom_smooth(aes(y = `European Union`), method = 'loess') +
  labs(x = 'Date',
       y = 'Issue attention',
       color = '',
       title = '') +
  scale_color_manual(values = colors) +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "bottom",
        legend.text = element_text(size = 15),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.margin=margin(0,0,0,0))

# saving plot for Figure 1 as 'figure_1.png'
ggsave("~/Replication materials/Figures/figure_1.png", 
       plot = p_vax_attention, 
      dpi = 600, 
       width = 8, 
       height = 7)


#### Figure 2: Interviews per day before and after the start date of the UK vaccine rollout ####

# deleting missing observations on treatment variable
data_filtered <- data |> 
  filter(is.na(treat2) == F) 

# building plot for Figure 2
p_fieldwork <- ggplot(data_filtered, aes(x = day)) + 
  geom_bar() +
  scale_x_continuous(limits = c(0, 31),
                     breaks = seq(0, 31, 5)) +
  geom_vline(xintercept = 17,
             color =  "red",
             linewidth = 0.5,
             linetype = "dashed") +
  labs(y = "Participants",
       x = "Day") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "bottom",
        legend.text = element_text(size = 15),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.margin=margin(0,0,0,0))

# Saving Figure 2 as 'figure_2.png'
ggsave("~/Replication materials/Figures/figure_2.png", 
       plot = p_fieldwork, 
       dpi = 600, 
       width = 8, 
       height = 7)


#### Table 1: Balance measures for covariate distributions ####

data_matching <- data |> 
  # Filter the dataset 'dat_2' to select the relevant variables of interest for analysis
  select(c("treat2", "country", "age", "gender", "eduy", "lr_self", "subjective_class", "living_area", "local_interest", "national_interest", "european_interest", "employment_r", "unemployment")) %>%
  # Remove any rows with missing values (NA) from the dataset
  na.omit() 

# Perform propensity score matching using the 'matchit' function
## 'treat2' is the treatment variable, and we match on the other covariates
m.out <- MatchIt::matchit(treat2 ~ age + gender + eduy + lr_self + subjective_class + living_area + local_interest + national_interest + european_interest + employment_r + unemployment,
                          data = data_matching) 

# Assess balance of covariates post-matching using the 'bal.tab' function
## The threshold for acceptable balance on the standardized mean difference is set at 0.1
balance <- bal.tab(m.out, thresholds = c(m = .1),
                   un = TRUE) 

# Convert the balance table into a data frame for easier manipulation and presentation
balancetable <- as.data.frame(balance$Balance) |> 
  rownames_to_column() |> # Move rownames into a new 'variable' column
  mutate(variable = rowname)

# Generate a summary of the matching output for both treated and control groups
summary <- summary(m.out,
                   un = TRUE)

# Convert the summary table into a data frame for easier manipulation and presentation
summarytable <- as.data.frame(summary$sum.all) |> 
  rownames_to_column() |> 
  mutate(variable = rowname)

# Combine relevant columns from both the balance table and the summary table
## This creates a final table comparing the balance metrics between treated and control groups
balance_test <- data.frame(Variable = balancetable$variable,
                           `Means treated` = summarytable$`Means Treated`,
                           `Means control` = summarytable$`Means Control`,
                           `SD mean difference` = summarytable$`Std. Mean Diff.`,
                           `Means threshold` = balancetable$M.Threshold)

# Display the final balance test results in a simple table format using 'kable'
t_balance <- kbl(balance_test,
                 digits = 4,
                 row.names = FALSE,
                 escape = FALSE)


save_kable(
  t_balance,
  "~/Replication materials/Tables/table_1.png",
  bs_theme = "journal",
  density = 600)


#### Figure 3: Effects of UK Vaccination on EU attitudes ####
## Building main models for DVs: EU support, EU beneficial, Positive image, EU right direction, EU democracy, Health decisions, EU crisis coordination, EU future crisis coordination
model_eu_support_fe <- lm_robust(eu_support ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                         data = data,
                                         fixed_effects = country,
                                         se_type = "stata")
model_eu_beneficial_fe <- lm_robust(eu_benefits ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                            data = data,
                                            fixed_effects = country,
                                            se_type = "stata")
model_eu_positive_fe <- lm_robust(eu_positive ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                          data = data,
                                          fixed_effects = country,
                                          se_type = "stata")
model_eu_right_direction_fe <- lm_robust(eu_right_direction ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                 data = data,
                                                 fixed_effects = country,
                                                 se_type = "stata")
model_eu_democracy_fe <- lm_robust(eu_democracy ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                           data = data,
                                           fixed_effects = country,
                                           se_type = "stata")
model_eu_health_issues_fe <- lm_robust(eu_health_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                               data = data,
                                               fixed_effects = country,
                                               se_type = "stata")
model_eu_future_coordination_fe <- lm_robust(eu_coord ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                     data = data,
                                                     fixed_effects = country,
                                                     se_type = "stata")
model_eu_coordination_fe <- lm_robust(eu_help ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                              data = data,
                                              fixed_effects = country,
                                              se_type = "stata")

## Creating a list of all main models
models_fe <- list("EU support" = model_eu_support_fe,
                  "EU beneficial" = model_eu_beneficial_fe,
                  "Positive image" = model_eu_positive_fe,
                  "EU right direction" = model_eu_right_direction_fe,
                  "EU democracy" = model_eu_democracy_fe,
                  "Health decisions" = model_eu_health_issues_fe,
                  "EU crisis coordination" = model_eu_coordination_fe,
                  "EU future crisis coordination" = model_eu_future_coordination_fe)

## Plotting main models 
p_main_models <- plot_summs(
  models_fe,
  inner_ci_level = .9,
  colors = brewer.pal(8, "Dark2"),
  coefs = c("Treatment" = "treat2")) +
  labs(y = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 13),
        legend.position = "bottom")

## Saving Figure 3 as 'figure_3.png'
ggsave("~/Replication materials/Figures/figure_3.png", 
       plot = p_main_models, 
       dpi = 600, 
       width = 9, 
       height = 7)

#### Figure 4: Treatment effects on attitudes towards European decision-making ####
## Building main models for DVs: EU health issues, social security, education, equality, climate, jobs, digitalisation, working conditions
model_eu_health_issues <- lm_robust(eu_health_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                              data = data,
                                              fixed_effects = country,
                                              se_type = "stata")
model_socialsecurity <- lm_robust(eu_ss_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                            data = data,
                                            fixed_effects = country,
                                            se_type = "stata")
model_education <- lm_robust(eu_education_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                       data = data,
                                       fixed_effects = country,
                                       se_type = "stata")
model_equality <- lm_robust(eu_equality_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                      data = data,
                                      fixed_effects = country,
                                      se_type = "stata")
model_climate <- lm_robust(eu_climate_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                     data = data,
                                     fixed_effects = country,
                                     se_type = "stata")
model_jobs <- lm_robust(eu_jobs_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                  data = data,
                                  fixed_effects = country,
                                  se_type = "stata")
model_digitalisation <- lm_robust(eu_digitalisation_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                            data = data,
                                            fixed_effects = country,
                                            se_type = "stata")
model_workingconditions <- lm_robust(eu_workingconditions_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                               data = data,
                                               fixed_effects = country,
                                               se_type = "stata")

## Creating a list of all specific variables models
models_specific <- list("Health" = model_eu_health_issues,
                                  "Social security" = model_socialsecurity,
                                  "Education" = model_education,
                                  "Equality" = model_equality,
                                  "Climate" = model_climate,
                                  "Working conditions" = model_workingconditions,
                                  "Jobs" = model_jobs,
                                  "Digitalisation" = model_digitalisation)

## Plotting specific variables models
p_models_specific <- plot_summs(
  models_specific,
  inner_ci_level = .9,
  title = 'Specific variables - portugal data FE',
  colors = brewer.pal(8, "Dark2"),
  coefs = c("Treatment" = "treat2")) +
  labs(y = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 13),
        legend.position = "bottom")

## Saving Figure 4 as 'figure_4.png'
ggsave("~/Replication materials/Figures/figure_4.png", 
       plot = p_models_specific, 
       dpi = 600, 
       width = 9, 
       height = 7)
